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The Magnet osp her ic MultiScale (MMS) mission employs a formation of spinning space- 
craft with several flexible appendages and thruster-based control. To understand the com- 
plex dynamic interaction of thruster actuation, appendage motion, and spin dynamics, 
each spacecraft is modeled as a tree of rigid bodies connected by spherical or gimballed 
joints. The method presented facilitates assembling by inspection the exact, nonlinear dy- 
namic equations of motion for a multibody spacecraft suitable for solution by numerical 
integration. The building block equations are derived by applying Newton’s and Euler’s 
equations of motion to an “element” consisting of two bodies and one joint (spherical and 
gimballed joints are considered separately). Patterns in the “mass” and “force” matrices 
guide assembly by inspection of a general N-body tree-topology system. Straightforward 
linear algebra operations are employed to eliminate extraneous constraint equations, re- 
sulting in a minimum-dimension system of equations to solve. This method thus combines 
a straightforward, easily-extendable, easily- mechanized formulation with an efficient com- 
puter implementation. 


I. Introduction 

T his paper presents the derivation of the exact nonlinear dynamic equations of motion for a multibody 
spacecraft with spherical or gimballed joints, suitable for computer solution. The formulation presented 
in this paper is motivated by the Magnetospheric MultiScale (MMS) mission, currently in mission definition 
at Goddard Space Flight Center. MMS will employ a formation of four spacecraft in a highly elliptical orbit 
to study the Earth’s magnetosphere. Each spacecraft is composed of a cylindrical central body, four “rigid” 
instrument booms, and four 50-meter wire booms. The spacecraft are spin-stabilized, and require periodic 
thruster operation for both formation maintenance and large orbital maneuvers. To study the interaction 
between thruster forces, spin dynamics, and boom flexibility, a simulation tool with sufficient fidelity is 
essential. To keep the model conceptually simple, we choose to model the booms (both rigid and wire) as 
chains of rigid bodies, with elastic and damping properties concentrated at the joints. To confirm the validity 
of this construction, we compare the simulated motion of one-, two-, and three-segment boom models; in 
the limit, the results should converge to that of the continuous structure. Seeking a dynamical formulation 
which allows easy adaptation from one-segment booms to two- and three-segment booms leads us to the 
general method presented here. 

The goal of this paper is to develop and describe algorithms in sufficient detail that the reader may, given 
properties of a multi-body system, assemble by inspection exact, nonlinear dynamical equations of motion 
suitable for numerical integration by computer. While the literature is rich with multi-degree-of-freedom 
dynamic formulations, most have theoretical or practical limitations which compromise their utility for 
problems like the MMS spacecraft. Finite-element models allow straightforward modelling of many modes; 
the method here presented borrows from finite-element methods in the assembly of large systems from small 
building blocks. Finite-element models, however, often linearize the equations of motion early in the process, 
and are consequently best suited to small perturbations from some nominal motion. Nonlinear coupling 
may be important to the MMS dynamics, and so a method was sought to preserve it. Hand derivation of 
nonlinear equations of motion using, for example, Kane’s or Lagrange’s methods, rapidly becomes laborious 
as the number of bodies grows. Computer aids help (Autolev, for example, even outputs source code), but 
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Figure 1. MMS Spacecraft, Scale and Closeup 


each model change (say from a one-segment boom to a three-segment boom) requires a new dynamic model 
to be linked to the simulation. The formulation presented in this paper allows the model information to be 
provided without recompilation of the simulation code. This point is exemplified by the simulation results 
presented, in which multiple spacecraft, with distinct dynamical models, are simulated concurrently using 
the same software functions. 

II. Equations of Motion for Two Bodies and One Spherical Joint 



Figure 2. Two Bodies Connected by a Spherical Joint 


Translational motion is governed by Newton’s second law, which we may write: 

mi) = F 


The rotational motion of a rigid body is governed by Euler’s equation: ’ 

Id) = T - u) x H 

where the overdot operator, ( ), denotes differentiaton with respect to time in a ‘local” frame; time derivatives 
of vectors in a frame A are related to derivatives in an inertial frame N by 

N -^-(v) = A -j-(u) + N lj a x v = i) -\- N u) A x v 
at at 
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For translational quantities, the most convenient “local” frame is the N frame itself. 

Consider the simplest possible multibody case, depicted in figure 2: an “inner” body, Bi , and an “outer” 
body, B 0 , connected by a spherical joint, G. We seek the equations of motion for attitude and translation 
of the two bodies with respect to an inertial reference frame, N. Considering each body separately, and 
explicitly separating the joint constraint force and torque from all other applied forces and torques, we may 
write Euler’s and Newton’s equations of motion as: 


li&i 

— Ti — LOi X Hi + Tq + Ti X Fq 

a) 

lo 

= T 0 — u 0 x H 0 — Tq — r Q x Fq 

( 2 ) 

rriiVi 

= Fi + Fg 

( 3 ) 

m 0 v 0 

$ 

1 

II 

( 4 ) 


where the symbols used for each body are listed in table 1, and where Fq and Tq are the constraint force 


Table 1. Symbols Used in Newton-Euler Equations 


Symbol 

Description 

LOi^LOo 

Angular velocity 

Vi , Vo 

Velocity of mass center 


Central moments of inertia 

mi,m 0 

Mass 

n,r 0 

Vector from mass center to joint 

H h H 0 

Central angular momentum 

Fi,F Q 

Resultant external force 

Ti, T 0 

Resultant external torque 


and torque at the joint. The local frames for differentiation of uji and u 0 are the frames fixed in Bi and B 0 , 
respectively; for both Vi and v Q , the N frame is the local frame. For a spherical joint, the constraint torque 
represents the material elasticity and damping lumped in the joint, as well as any active control torques; the 
constitutive relations to compute Tq are assumed to be provided by the user. The constraint force, Fb, is 
an unknown to be eliminated. We arbitrarily choose the signs of Fb and Tq to act in a positive sense on 
Bi , and in a negative sense on B 0 . In the absence of momentum storage devices, H — Iu. For completeness, 
the addition of embedded momentum wheels is addressed in the Appendix. 

The joint constraint is introduced by equating two expressions for the velocity of G in N: 


vg = Vi + ui x n = v 0 + u 0 x r Q 

Differentiating with respect to time (in N) yields the constraint equation: 

ii+WjXri+WiX (ui x n) =i' 0 +w 0 xr 0 + o; 0 x (u 0 x r Q ) ( 5 ) 

Equations (l)-( 5 ) are five vector equations in the five vector unknowns a;*, u Q , i v Q , and Fb. 

For compactness, we introduce two notational devices. The over-tilde operator, (), denotes a 3 -by -3 
skew-symmetric matrix so that xy = x xy. Similarly, the overbar operator, ( ), denotes a symmetric matrix 
such that xy = x x (x x y). If the scalar components of the vector x are (xi, x 2 , £3)5 then 



0 

-X3 

X2 


rr 2 

x 2 x 3 

X\X2 

X1X3 

x = 

£3 

0 

—X\ 

, x = 

X 2 Xi 

rp 2 ^,2 

x 3 x 1 

X2X3 


. ~ x 2 

Xi 

0 


X3X1 

X3X2 

rp 2 rp 2 

x 1 x 2 


Note that —x = x T . 
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Now we write equations (l)-(5) in matrix form: 


h 

0 

0 

0 

-n " 


• N 

Ui 


Ti - CjiHi + T g 'j 

0 

Io 

0 

0 

To 


UJo 


1 

0 

1 

0 

0 

mil 

0 

-l 

< 

Vi 

> = < 

F < 

0 

0 

0 

m 0 1 

l 


Vo 


Fa 

~T 

-r\ 

~T 

r o 

-1 

1 

0 


:f g , 


„ - w 0 r 0 


where 1 denotes the 3-by-3 identity matrix. 


Table 2. Reference Frames for Expression of Components of Dynamic Entities 


Items Expressed in Bi 

Items Expressed in B 0 

Items Expressed in TV 

h 

Io 

Vi, Vi 

n, n 

r 0 , r Q 

Vot V Q 

UJi, LJi 

tJo) 

Fi 

Ti 

T 0 

F 0 

Hi 

H 0 

f g 

t g 




( 6 ) 


These equations of motion are basis-free. To solve them numerically, all items in an equation must be 
expressed in a common reference frame. Table 2 lists dynamic entities of interest and the frames in which 
they are naturally expressed. Suppose, for example, that a thruster is mounted on Bi. The force it exerts 
on Bi is expressed in the TV frame for addition into F*, while the torque it exerts on Bi is expressed in Bi , 
as it must be added into T*. It proves most convenient to express the first row of equation (6) in J5*, the 
second in J5 0 , and the remaining three in TV. The choice to express Tq in Bi rather than in B 0 is arbitrary. 

The transformation between reference frames X and Y is represented by the direction cosine matrix X C Y . 
If a vector v is expressed in the Y frame as Y v , its components in the X frame are computed by 

Note that y C x = ( X C Y ) T = ( X C Y )~\ 

Introducing bases, we obtain 


Ii 

0 

0 

0 

-fi V N ' 


f • > 

Ui 


' T 

- QiHi + T G ' 

0 

Io 

0 

0 

fo°C N 


U 0 


T 0 - 

- Cj 0 H 0 - °C i T G 

0 

0 

TO jl 

0 

-1 

< 

Vi 

> = < 


Fi 

0 

0 

0 

TO 0 1 

1 


Vo 



F 0 

. - N C l fJ 

N C°r T 0 

-1 

1 

0 


< Fg > 


k N C { 

UiTi - N C 0 w 0 r 0 . 


This equation is the building block for the multibody case. When assembled, the system equations of motion 
will be of order 67V + 3 (TV — 1), representing three rotational degrees of freedom per body, three translational 
degrees of freedom per body, and three constraints for each joint. We will see in section IV how to reduce 
the number of equations to a minimal set. 

III. Generalization to N Bodies 

This section develops the procedure for assembling the system equations of motion using equation (7) as 
the building block. A tree topology of TV bodies has TV — 1 joints. While a body may be associated with 
multiple joints, each joint has exactly one “inner” body and one “outer” body. Using this as an organizing 
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principle allows the assembly of the system equations of motion by inspection. First, we write: 






f 

LOi 

h 

0 

p 


UN 

In 

mil 


< 

Vi 

: 

0 


j 


V N 


m^v 1 



Fgi 

P T 

j t 

0 


\ 





< f g{n- i) , 


where 


Ti — Cb\H\ + y ^Tc(ji) 


Tn — &nHn + ^2 Tc(jN ) 


Fn 

- { N C°uJ 0 r 0 )Gi 


[ ( N C l UJiri)G(N- 1 ) — ( N C O L0 o ro)G{N-l) 

( 8 ) 


T G{jk ) = \ 


Tg 

- °C i T c 
0 


if body k is Bi for joint Gj 
if body k is B 0 for joint Gj 
otherwise 


( 9 ) 


> 


The Tc(jk) terms are assembled by a simple procedure. For each joint, apply the (user-supplied) joint torque 
to the joint’s inner body, and apply its equal and opposite (properly transformed) to the joint’s outer body. 

The submatrices P and J share a common structure. Each may be considered as an array of 3 x 3 blocks, 
arranged in N “3-rows” and N — 1 “3-columns”. Each such 3-row corresponds to a body, and each 3-column 
corresponds to a joint. Thus, each 3-column contains exactly two nonzero blocks, corresponding to the inner 
and outer bodies associated with that joint. For each joint Gj, 


P ij = -f ij i C N , P 


oj 


Jij — 1 5 


' OJ 

Joj — 1 


°c N 


( 10 ) 


As an illustrative example, consider the system depicted in figure 3. We construct a tree connectivity 
table, assigning the inner and outer body for each joint. See table 3. 



Figure 3. Four-body Example Problem 


The submatrices P and J follow the structure of the tree connectivity table: 



p 

1 

0 

0 


’ -1 

0 

0 

p = 

f 21 2 C N 

0 

-f 22 2 C N 

T32 3 C N 

-r 23 2 C N 

0 

, J = 

1 

0 

-1 

1 

-1 

0 


0 

0 

f 43 4 C" . 


0 

0 

1 
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Table 3. Example Tree Connectivity Table 



G 1 

g 2 

g 3 

B 1 

l 



b 2 

0 

i 

i 

B 3 


0 


B 4 



0 


The right-hand side of equation (8) is written explicitly as 

T\ — w\H\ + T G i 

T2 — OJ2H2 — 2 C 1 Tgi + Tq 2 + Tg 3 
T 3 -Cj 3 H 3 - 3 C 2 Tg2 
T 4 -u a H 4 - a C 2 T G 3 
Fi 

< F 2 > 

F 3 
F a 

N C 1 ZJirn - N C 2 oJ 2 r2 1 
N C 2 uJ 2 r 22 — N C 3 uJ 3 r 32 
N C 2 uJ 2 r 23 - N C 4 ZJaU 3 


IV. Reduction of the State Vector 


The equations of motion described in the previous section have a clear structure deriving from their 
origin in Newton and Euler’s equations. This structure facilitates mechanization; the various submatrices 
may be readily assembled by simple procedures. This formulation suffers, however, from rapid growth of 
computational effort with the size of the model. The number of equations, 6N + 3(N — 1), grows three times 
as fast as the number of degrees of freedom, 3N + 3. To reduce the computational load, we assume that the 
joint forces, are not of interest, and seek a reduced, minimum-dimension system of 3N + 3 equations to 
solve. 

We begin by defining a reduced state vector x, which has all the desired degrees of freedom, and two 
other auxiliary vectors y and /: 


x = 


Ui 1 


W 2 


f N 

V2 


( 

Fg 1 


> , y=< 


►, f=< 

: > 

UN 


VN 

\ > 


, F G(N- 1) , 


Vi ) 


We then repartition equation (8) thus: 


A 

0 

R " 


f X 1 


f T ' 

0 

M 

U 

< 

1 . 1 
1 y 1 

► = < 


R T 

U T 

0 _ 


l/J 


[c J 


where the system matrix partitions are correspondingly: 


(ii) 
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T= < 



" II 


|- -i 




7712 1 

A = 

Iat 

, M = 

771 AT 1 


mil 




t 1 -q 1 h 1 + Y,Tgud ' 





3 


f > 

f 2 


( N C'uiri)Gi ~ ( N C°u) 0 r 0 ) G i 

Tn — &nHn + ^ Tg(jn) 

i 

> , T — < 


► , c= < 

: > 
( N C l uJiri)o(N-i ) ~ { n C°U 0 '>’o)g(n- i) _ 


Fi 


R is formed by adjoining the first 3-row of J onto P, and U is the portion of J that remains. In our 
four-body example, 


R = 


-r 11 l C N 0 0 

r 2 i 2 C N -v 2 2 2 C n -r 2 s 2 C N 
0 r 32 3 C N 0 

0 0 f 43 4 C N 

-10 0 


U = 


1 

0 

0 


-1 

1 

0 


-1 

0 

1 


We now perform row operations to decouple x from y and /. Introducing the unknown coefficient matrices 
a and /?, we add the three 3-rows of equation (11): 


A 0 R 


+ a 


0 M [/]+/?[ R t U t 0 ]] < y 


> = T + aT + PC 


or 


/ \ 

x 


[ A + /3R T aM + /3U T R + aU ] < y 

l f J 

To decouple y and /, we solve for a and (3 so that 

aM + (3U t = 0 
R odJ — 0 


> — T + olT + PC 


( 12 ) 


from which we find 


a = —RU~ l 
(3 = RU~ l MU- T 


and the desired form of the equations of motion is 

[A + (3R t ) x = T + aT + f3C (13) 

This is the system of equations to be solved at each step of the numerical integration. Note that A, M, 
and U, (and thus U~ l and U~ 1 MU~ T ) are constant matrices, and so need only be computed once for a 
particular system. 
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V. Equations of Motion of Two Bodies with One Gimballed Joint 


To this point, the dynamical formulation assumed spherical joints, possessing three rotational degrees of 
freedom with no preferred rotation axis. In spacecraft dynamic modelling, applications commonly arise where 
it is desirable to model a joint having one or two degrees of freedom as a gimballed joint; solar array drives 
and antenna pointing mechanisms are prime examples. So, in the following sections, we develop equations of 
motion for a system of rigid bodies linked by gimballed joints, building on the framework already developed 
for spherical joints. 

The rotation of a gimbal joint is described by an Euler sequence of gimbal angles. Considering again 
the system described by equations (1)— (5), we introduce the gimbal angles, 0, and their associated rates of 
change, u = 0, so that the angular velocity of the outer body may be written as 

u 0 = uJi + Tu (14) 

where T is a 3-by-3 matrix of nonlinear functions of the gimbal angles, depending on the given Euler 
sequence and whether the rotations are about “space” or “body” axes; see Kane, Likens, and Levinson 1 
for a complete discussion of these terms. For example, a Body- three 1-2-3 Euler rotation through 0 1, 0 2 , #3 
has the relationship between 1 uj° and u : 


i, ,o | 

V 1 


COS 02 COS 03 

sin 03 

o ’ 


r ui 

l(jJ 2 | 

► = 

— cos 02 sin 03 

COS 03 

0 

< 

| U2 > 

luJ 3 J 


sin 02 

0 

1 




Writing it compactly as 1 uj° = Tu defines T for this Euler rotation sequence. Differentiating (14), 

uj 0 — Cji -\- r ii -|- r u 


where f is, like T, a nonlinear function of 0 and u. Substituting into (5), some manipulation yields 


Vo = Vi + (r Q - Vi) LJi + To (tu + tu'j + ujiTi - UJ 0 T 0 

and the equations of motion for the two-body system are 



0 

0 

-n 


LOi 


Ti — UiHi + Tq 

Io 

io r 

0 

r Q 

< 

it 

> — > 

T 0 — u> 0 H 0 — Tg — I 0 tu 

0 

0 

mil 

-1 


Vi 

r — > 

Fi 

. m 0 (r Q - fi) 

m 0 f o r 

m Q 1 

1 


< f g j 


k F 0 ~ m 0 f 0 tu + m Q (ujiTi - uJ 0 r 0 ) , 


(15) 


Comparing with equation (6), we see some familiar patterns, but not enough to extrapolate with confidence 
to the N-body problem. For further illumination, we turn to the four-body example of figure 3. The equations 
of motion for this example may be written 


I 0 P 

n fi j 


u 1 

iii 
U 2 

l ^3 ) 


Vi ! 

F g ) 



(16) 


with P and J as defined in equation (10) except that we have not yet introduced coordinate bases, and 


' h 

0 

0 

0 


mil 

h 

I 2 T 1 

0 

0 

II — 

m 2 l 

h 

I 3 T 1 

I 3 T 2 

0 

5 t 1 

m 3 l 

. Ii 

hTi 

0 

I 4 T 3 . 


7714! 


8 of 13 


American Institute of Aeronautics and Astronautics 



n = 


0 

0 

0 

0 

m 2 (r 2 1 - rn) 

m 2 r 2 iTi 

0 

0 

m 3 (f 32 - r 22 + f 2i - rn) 

m 3 (f 32 -f 22 +f 2 i)]?i 

m 3 f 32 r 2 

0 

m 4 (r 43 - f 23 + r 2i - fn) 

mi (f 43 - r 23 +f 2 i)Ti 

0 

m 4 f 43 r 3 _ 


T 


Ti — £j\H\ + Tgi 

T 2 — 0J2H2 — Tgi 4 * Tg2 + Tg 3 — 7 2 f i u i 
Ts — U3H3 — Tg2 ~ h (fitii 4 * f 2U2) 

T4 — — Tg 3 — ^4 (r 1^1 4 * T3?i3j 

Fl 

F 2 - m 2 f 2 if 1U1 4 - m 2 (oJirn - u; 2 r 2 i) 

F 3 - m 3 (f 32 - f 22 4 - r 2 i) 14^1 — m 3 f3 2 f 2 u 2 + m 3 (oJirn - aJ 2 r 2 i) + ra 3 (o; 2 r 22 - U3V32) 
F4 — m\ (7*43 — f 2 3 + r 2 i) f \U\ — m^r^tsus + 7714 (aJi rn — uJ 2 r 2 i) + 7714 (uJ 2 r 2 3 — ^4^43) 


While some of the symmetries of the spherical joint system are lost, new patterns emerge. The I and II 
matrices have a common structure stemming from the composition of angular velocities by concatenation of 
gimbal rates with uj For a system of N bodies, the matrix I is 3N x 3iV, each 3-row j corresponding to a 
body, and each 3-column k corresponding to a point. The first 3-column corresponds to the mass center of 
Bi, since it multiplies Cj\. Each succeeding 3-column corresponds to a joint. Each 3-row may be assembled 
by inspection, by traversing the tree from B\ to Bj , placing an IjTk-i term for each joint in the path. 
Similarly, the elements of the II matrix are parallel-axis (“ mr 2 ”) terms, the r vector being the vector from 
the mass center of Bj to point k. The novel terms in T and T may be likewise constructed; compare /, 
II, T, and T with the path traversals tabulated in table 4. These patterns will lead us to codify rules for 
assembly of equations of motion for a general system. First we must introduce coordinate bases. 


Table 4. Example Tree Path Table 



G 1 

G 2 

G3 

B 1 




b 2 

X 



b 3 

X 

X 


Ba 

X 


X 


Assigning bases proceeds as with the spherical joints, except that Tit and f u terms are naturally expressed 
in the inner body frame for “space” sequences, but in the outer body for tc body” sequences. Subsequent 
discussion assumes a body sequence. 

Considering again the two-body case, we observe that it is most convenient to express the first equation 
in Bi , the second in B 0 , and the remaining two in the N frame. Inserting the appropriate direction cosine 
matrices, equation (15) becomes 


It 

0 0 

1 

e*. 

0 

1 


( 

Io 

Io r 0 

r 0 °C N 

< 

u 

0 

0 mil 

-1 


Vi 

. m 0 ( N C°f 0 - N C i f i ) 

m 0 N C°f 0 r m G 1 

1 


. Fg 


Ti — uJiHi + T g 
T q — Q 0 H 0 — o&Tg-IoTu 
Fi 

< F 0 — m 0 N C°r 0 tu + m 0 ( N C i u i r i - N C°cJ 0 r 0 ) j 


(17) 


The equations of motion for the assembled system are still equation (16). For notational convenience, we 
consider the mass center of B\ to be the zeroth joint, Go, so that To = 1 and ( 7^)0 = 0. Then the system 
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submatrices may be written: 


Ujk = rrij 


Ijk — 


( N c°f 0 ) 


k - 1 


IjTk-i if Gk-i is in path from B\ 

0 otherwise 

+ ^2( N C°f 0 - N C i fi i C°) 

1 



to Bj 


3,k = 

l : All 


j,k = l,...,N 
i 

G in path between Gk-i and Bj 


Hj = rrijl, 


j = 1, . . . , N 


P and J are defined as in equation (10). 


j,k = l,...,N 

Tj = Tj — CjjHj + ^2 ^G(jfc) — Ij j C l t[Ui Tg defined in equation (9) 

fe 1 l : All G in path from B\ to Bj 



( N c°n) k + Y,C' c ° f ° 

l 


N C i f i i C°) l 


t k u k + J2[ N C 0 Lj 0 ro 

k 


N 


C l u)iri 


) 


j = l,...,N 

k : All G in path from B\ to Bj 
l : All G in path between Gk and Bj 


VI. Elimination of Constraint Forces and Constrained Gimbal Degrees of 

Freedom 


Recall from section IV that the spherical joint formulation led to a system of dimension 9 N — 3, which 
we then reduced to a 3N -f 3 system by row operations. The gimballed joint formulation starts from a more 
advantageous position, in that the redundant velocity and angular velocity states have already been reduced 
by combination into the gimbal angular rates, it, resulting in a system of dimension 6N. To arrive at a 3Af + 3 
system, we eliminate the 3 N — 3 joint forces, Fq, by row operations. Further, the motivation for pursuing 
the gimballed formulation is to eliminate undesired degrees of freedom, such as the two constrained axes of 
a single-axis solar array drive mechanism, from the final equations of motion. With this formulation, these 
degrees of freedom are easily eliminated. 

First we eliminate the joint constraint forces by defining the state vector as 


Ui 

u 

Vl 


and partitioning the assembled equations of motion as 


A R 
S U 


U) 


T 

P 


The subsystem matrix A is formed by adjoining the first 3-rows of II and /( to the 7 0 of equation (16): 


A = 


I 0 
III Mi 


and S comprises the remainders of II and /x: 


S — 112,... ,N M2,...,AT 
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As for the spherical-joint case, the submatrix R is formed by adjoining the first 3-row of J onto P, and U is 
the portion of J that remains. Again using an unknown coefficient matrix, a , so that 


(A + aS) x = T + olT 

R -I- olU — 0 


we solve for ct, 


a = -RU- 1 


Substituting, we obtain the reduced equations of motion, 


(A - RU~ 1 S) x = T - RU' 1 ? 


( 18 ) 


Finally, we consider the gimbal degrees of freedom which are to be eliminated due to their being con- 
strained or negligible. In the spherical joint case, these degrees of freedom were complicated combinations 
of the states. With the dynamics reformulated in terms of gimbal states, eliminating undesired degrees of 
freedom is as straightforward as striking out rows and columns of the system matrices. To make this clear, 
expand the state vector into its scalar components, 

x — { Un U >12 to 13 Un U \2 ^13 • • • UNI UN2 } 

where Ujk is associated with the kth degree of freedom of the j th joint. Suppose, as an example, that joint 
1 is a two-axis gimbal, so that u\s = uis = 0. Looking at the system of equation (18), the sixth column of 
(A — RU~ 1 S) multiplies zero, and the sixth row becomes a constraint equation on the reaction torque at 
the gimbal. Eliminating the sixth row and column of the system removes the undesired degree of freedom, 
reducing the system order by one. In the implementation, this state reduction may be accomplished during 
the assembly of A,P, 5, and T. 


VII, Some Simulation Results 

As a demonstration of the methods developed in this paper, some early MMS simulation results are 
presented. As has been discussed, the interaction of spin, boom flexibility, and thruster activity provides 
some interesting problems. A thorough discussion of these is beyond the scope of this paper; our purpose 
here is simply to show that the numerical methods have been successfully implemented. 

Three spacecraft models are implemented. Each consists of a central body with four “wire” appendages. 
For the first spacecraft, each appendage is modeled as a single rigid body. The second and third spacecraft 
model each appendage as two and three connected rigid bodies, respectively. Spherical joints are used, 
with joint torques consisting of light damping provided through a user-supplied function. The numerical 
integration of the dynamics is accomplished using the same software functions for all three models. The body 
and joint information that defines each model is input from text files and stored in assigned data structures, 
and the dynamic system matrices are assembled as needed at each simulation timestep. All three spacecraft 
models are simulated concurrently. 

The booms are given small initial deflections to excite both in-plane and out-of-plane oscillations. Figure 
4 shows the motion of one boom tip, for all three models. All models show the interaction between in-plane 
and out-of-plane boom dynamics, coupled through spacecraft spin. The second and third modes of the boom 
are most evident in the initial transient of the in-plane tip motion of the two- and three-segment models. 
These higher-frequency modes are attenuated by the damping torques in the joints. The model comparison 
lends confidence both that the important dynamic features have been captured in the model, and that the 
mathematical models have been correctly implemented in software. 

VIII. Conclusion 

Beginning from Newton’s and Euler’s equations of motion for a single rigid body, we have presented a 
method for constructing the dynamic equations of motion of a general system of rigid bodies connected by 
spherical or gimbaflled joints in a tree topology. Two guiding ideas make this a tractable problem. First, each 
joint in a tree topology has exactly one ‘inner” body and one “outer” body. This allows a system of arbitrary 
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Figure 4. Wire Boom Tip Motion 

complexity to be assembled from the equations of motion of a two-body system. The assembly procedure 
is easy to implement in software, and relatively easy to check by inspection of the various arrays as they 
are populated. The second idea is that elimination of constraint equations may be considered as a problem 
in linear algebra to be solved numerically at each timestep, rather than a dynamics problem to be solved 
during formulation. For simple systems, it may be practical to symbolically eliminate the 3 N — 3 constraint 
equations, yielding some analytical insight. This rapidly becomes impractical, however, as the system grows 
in complexity. The assembly method presented here allows straightforward verification that the matrices 
A , M, etc. have been properly loaded into memory, as the patterns and symmetries in the building blocks 
are not obscured. Elimination of the constraint equations is then a purely mechanical process of matrix 
inversions and multiplication. 

Appendix: Auxiliary Equations of Motion 

Kinematic Equations of Motion 

The main body of this paper discusses the dynamic equations of motion, which are solved to find the dynamic 
state variables, lj and v. For completeness, we present without derivation the additional equations of motion 
associated with the kinematic state variables: the quaternion, q , representing attitude, and the position 
vector, p. 

Defining the quaternion elements in terms of a simple rotation about an axis A through an angle 0: 

qi ( Aisin(0/2) 1 

< q2 > = i sin ( 0 / 2 ) i 

q 3 | A 3 sin(<9/2) | 

< q 4 ) [ cos(<9/2) J 

the time derivatives of the quaternion elements are functions of the quaternion and the angular velocity: 



The translational kinematic equations of motion are simply: 

p — v 

Equations of Motion for Embedded Momentum Wheels 

For a single rigid body B , Euler’s equation is written: 

Iuj = T - lj x H 
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where H — Iuj. We consider here the case of one or more momentum wheels embedded in a rigid body. 
Consider a reaction wheel, W , having rotational inertia J, scalar angular momentum h w , and whose axis, 
a, is fixed in B. Momentum is exchanged between B and W by applying equal and opposite interaction 
torques (typically specified by some control law). We denote this interaction torque as a scalar, T w m , we apply 
it (arbitrarily) in a positive sense to the wheel, and its opposite to B. The modified equations of motion are: 

Iuj — T — uj x (/cj + h w a) — T w a 

h , 'iu — '-L \jj J d * uj 

Multiple wheels act in superposition: 

n n 

IuJ — T UJ X (Jo; *4“ ^ ^tu/e^/c) ^ ^ T w kCLk 

k = 1 /c=l 

h"wk = T-'wk Jkdk * ^5 (k = 1 , 2 ,..., ? 7 <) 
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